library(ArchR)
library(rhdf5)
library(reticulate)
library(tidyverse)
#library(caret)
h5disableFileLocking()
# load ArchR project
proj <- loadArchRProject(path = "Save_Vignette_object_4/")
# for R to find macs2 we will use a conda environment where the latest version
# of macs2 is installed
#conda_list()[[1]][2] %>% use_condaenv(required = TRUE)
If we compute correlation on sparse features, this can lead to noise in the data. Therefore, first low-overlapping aggregates of single cells are created (see Cicero). Aggregates with >80% overlap with an other aggregate are removed.
Correlation in accessiblity between two peaks across many single cells. For example if peak A is accessible in a single cell is peak B also accessible? With this analysis we might learn if a certain enhancer peak is often co-accessible with a certain promoter peak in a particular cell type.
Keep in mind that cell-type specific peaks would be identified as co-accessible within this cell tpye. This strong correlation does not mean that there is an underlying co-accessibility!
Query hits and subject hits are the indices of the two peaks that are correlated. The correlation column contains the numeric correlation value of the accessibility between the two peaks.
# compute peak co-accessibility information and store it in the ArchR project
proj <- addCoAccessibility(proj,
reducedDims = "IterativeLSI",
k = 100,
corCutOff = 0.75,
# max. 80 % overlap between current group and all previous groups to permit
# current grup to be added to the group list during k-nearest neighbor calculation
overlapCutoff = 80,
dimsToUse = 1:30,
#number of k-nearest neighbor groupings to test for passing the supplied overlapCutoff
knnIteration = 500 )
# to retrieve the co-accessibility information we use getCoaccessibility()
# return dataframe if returnLoops = FALSE
coacc <- getCoAccessibility(
ArchRProj = proj,
# if the dimension has a correaltion to the sequencing depth
# that is greater than the cutoff it will be excluded
corCutOff = 0.5,
resolution = 1,
returnLoops = FALSE # if TRUE it returns a GRangesObject
)
coacc %>% head %>% knitr::kable()
| queryHits | subjectHits | seqnames | correlation | Variability1 | Variability2 | TStat | Pval | FDR | VarQuantile1 | VarQuantile2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 6 | 9 | chr1 | 0.5948088 | 0.0041591 | 0.0302828 | 16.51231 | 0 | 0 | 0.4991376 | 0.9302221 |
| 9 | 6 | chr1 | 0.5948088 | 0.0302828 | 0.0041591 | 16.51231 | 0 | 0 | 0.9302221 | 0.4991376 |
| 9 | 31 | chr1 | 0.6909559 | 0.0302828 | 0.0141892 | 21.32989 | 0 | 0 | 0.9302221 | 0.8024554 |
| 18 | 19 | chr1 | 0.5128650 | 0.0259368 | 0.0052909 | 13.33193 | 0 | 0 | 0.9060243 | 0.5682542 |
| 19 | 18 | chr1 | 0.5128650 | 0.0052909 | 0.0259368 | 13.33193 | 0 | 0 | 0.5682542 | 0.9060243 |
| 21 | 49 | chr1 | 0.5345478 | 0.0042255 | 0.0335266 | 14.11476 | 0 | 0 | 0.5038296 | 0.9449599 |
Metadata:
metadata(coacc)
## $peakSet
## GRanges object with 141025 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## Mono chr1 752506-753006 *
## B chr1 762688-763188 *
## GMP chr1 773653-774153 *
## B chr1 801007-801507 *
## B chr1 805041-805541 *
## ... ... ... ...
## Mono chrX 154840753-154841253 *
## NK chrX 154842390-154842890 *
## NK chrX 154862033-154862533 *
## Erythroid chrX 154912392-154912892 *
## GMP chrX 154997007-154997507 *
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
dim(coacc)
## [1] 111132 11
We can also return the co-accessibility information as a loop track. The resolution parameter sets the base-pair resolution of the loops. Resolution = 1 creates loops that connect the center of each peak.
Where is the information contained which peaks are correlated?
coacc <- getCoAccessibility(ArchRProj = proj,
corCutOff = 0.5,
resolution = 1,
returnLoops = TRUE)
coacc [[1]]
## GRanges object with 55566 ranges and 9 metadata columns:
## seqnames ranges strand | correlation Variability1
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 845703-856628 * | 0.594809 0.00415906
## [2] chr1 856628-940499 * | 0.690956 0.01418923
## [3] chr1 894695-895356 * | 0.512865 0.00529094
## [4] chr1 901579-999702 * | 0.534548 0.00422549
## [5] chr1 948868-1004761 * | 0.506462 0.02580794
## ... ... ... ... . ... ...
## [55562] chrX 153583564-153665944 * | 0.537871 0.000662602
## [55563] chrX 153584852-153597227 * | 0.566064 0.001244477
## [55564] chrX 153597227-153637613 * | 0.512202 0.063475852
## [55565] chrX 153686270-153744595 * | 0.513243 0.012864089
## [55566] chrX 153947827-154027599 * | 0.524475 0.013135018
## Variability2 TStat Pval FDR VarQuantile1
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 0.0302828 16.5123 3.62193e-49 2.50640e-47 0.499138
## [2] 0.0302828 21.3299 3.41000e-72 6.81189e-70 0.802455
## [3] 0.0259368 13.3319 6.84909e-35 2.23707e-33 0.568254
## [4] 0.0335266 14.1148 2.71023e-38 1.06988e-36 0.503830
## [5] 0.0189477 13.1076 6.23013e-34 1.92659e-32 0.905221
## ... ... ... ... ... ...
## [55562] 0.00116998 14.2381 7.75355e-39 3.14586e-37 0.0623195
## [55563] 0.06347585 15.3237 1.05928e-43 5.56776e-42 0.1644870
## [55564] 0.02888357 13.3086 8.62737e-35 2.80265e-33 0.9912758
## [55565] 0.02225960 13.3453 6.00370e-35 1.96741e-33 0.7838576
## [55566] 0.00519124 13.7465 1.10713e-36 3.99906e-35 0.7879772
## VarQuantile2 value
## <numeric> <numeric>
## [1] 0.930222 0.594809
## [2] 0.930222 0.690956
## [3] 0.906024 0.512865
## [4] 0.944960 0.534548
## [5] 0.853012 0.506462
## ... ... ...
## [55562] 0.151581 0.537871
## [55563] 0.991276 0.566064
## [55564] 0.923357 0.512202
## [55565] 0.880343 0.513243
## [55566] 0.562910 0.524475
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
If we decrease the resolution to 1000 this could help with overlplotting of co-accessibility interactions. The GRanges object will contain fewer total entries, 51,82 instaed of 54,002.
cA <- getCoAccessibility(
ArchRProj = proj,
corCutOff = 0.5,
resolution = 1000,
returnLoops = TRUE
)
cA[[1]]
## GRanges object with 53116 ranges and 9 metadata columns:
## seqnames ranges strand | correlation Variability1
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 845500-856500 * | 0.594809 0.00415906
## [2] chr1 856500-940500 * | 0.690956 0.01418923
## [3] chr1 894500-895500 * | 0.512865 0.00529094
## [4] chr1 901500-999500 * | 0.534548 0.00422549
## [5] chr1 948500-1004500 * | 0.506462 0.02580794
## ... ... ... ... . ... ...
## [53112] chrX 153583500-153665500 * | 0.537871 0.000662602
## [53113] chrX 153584500-153597500 * | 0.566064 0.001244477
## [53114] chrX 153597500-153637500 * | 0.512202 0.063475852
## [53115] chrX 153686500-153744500 * | 0.513243 0.012864089
## [53116] chrX 153947500-154027500 * | 0.524475 0.013135018
## Variability2 TStat Pval FDR VarQuantile1
## <numeric> <numeric> <numeric> <numeric> <numeric>
## [1] 0.0302828 16.5123 3.62193e-49 2.50640e-47 0.499138
## [2] 0.0302828 21.3299 3.41000e-72 6.81189e-70 0.802455
## [3] 0.0259368 13.3319 6.84909e-35 2.23707e-33 0.568254
## [4] 0.0335266 14.1148 2.71023e-38 1.06988e-36 0.503830
## [5] 0.0189477 13.1076 6.23013e-34 1.92659e-32 0.905221
## ... ... ... ... ... ...
## [53112] 0.00116998 14.2381 7.75355e-39 3.14586e-37 0.0623195
## [53113] 0.06347585 15.3237 1.05928e-43 5.56776e-42 0.1644870
## [53114] 0.02888357 13.3086 8.62737e-35 2.80265e-33 0.9912758
## [53115] 0.02225960 13.3453 6.00370e-35 1.96741e-33 0.7838576
## [53116] 0.00519124 13.7465 1.10713e-36 3.99906e-35 0.7879772
## VarQuantile2 value
## <numeric> <numeric>
## [1] 0.930222 0.594809
## [2] 0.930222 0.690956
## [3] 0.906024 0.512865
## [4] 0.944960 0.534548
## [5] 0.853012 0.506462
## ... ... ...
## [53112] 0.151581 0.537871
## [53113] 0.991276 0.566064
## [53114] 0.923357 0.512202
## [53115] 0.880343 0.513243
## [53116] 0.562910 0.524475
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
We can add the co-accessibility information as loop tracks to our browser track. For this we can use the loops parameter.
# marker genes we are interested in
markerGenes <- c(
"CD34", #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", #B-Cell Trajectory
"CD14", #Monocytes
"CD3D", "CD8A", "TBX21", "IL7R" #TCells
)
# plot
p <- plotBrowserTrack(
ArchRProj = proj,
groupBy = "Clusters2",
geneSymbol = markerGenes,
upstream = 50000,
downstream = 50000,
loops = getCoAccessibility(proj)
)
## GRanges object with 9 ranges and 2 metadata columns:
## seqnames ranges strand | gene_id symbol
## <Rle> <IRanges> <Rle> | <character> <character>
## [1] chr1 208059883-208084683 - | 947 CD34
## [2] chrX 48644982-48652717 + | 2623 GATA1
## [3] chr9 36838531-37034476 - | 5079 PAX5
## [4] chr11 60223282-60238225 + | 931 MS4A1
## [5] chr5 140011313-140013286 - | 929 CD14
## [6] chr11 118209789-118213459 - | 915 CD3D
## [7] chr2 87011728-87035519 - | 925 CD8A
## [8] chr17 45810610-45823485 + | 30009 TBX21
## [9] chr5 35856977-35879705 + | 3575 IL7R
## -------
## seqinfo: 24 sequences from hg19 genome
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
grid::grid.newpage()
grid::grid.draw(p$CD14)
What is the difference between co-accessibility and peak2gene linkage?
Co-accessibility uses only scATAC-seq data looking for correlations in accesssiblity.
Peak2Gene linkage in comparison also uses integrated scRNA-seq data to look for correlations between peak accessibility and gene expression.
Both approaces are solutions to a similar problem, but peak2gene linkage is assumed to provide more biologically relevant regulatory links.
All potential peak-to-gene linkages are identified and significant links (R > 0.45 and FDR < 0.1) are kept (Higher correlation and variance suggest that a link is not random noise.)
proj <- addPeak2GeneLinks(
ArchRProj = proj,
reducedDims = "IterativeLSI",
useMatrix = "GeneIntegrationMatrix",
dimsToUse = 1:30,
k = 100, # number of k-nearest neighbors to use for creating single cell goups for correlation analysis
maxDist = 250000, # maximum allowable distance in basepairs between two peaks to consider for co-accessiblity
scaleTo = 10^4 # total insertion coutns from a group of cell is summed across all peaks and normalized to toal depth provided by scaleTo
)
p2g <- getPeak2GeneLinks(
ArchRProj = proj,
corCutOff = 0.45,
resolution = 1,
returnLoops = FALSE,
varCutOffATAC = .25, # minimum variance quantile of ATAC peak accessibility
varCutOffRNA = .25, # minimum variance quantile of the RNA gene expression when selecting links
FDRCutOff = 1e-04 # maximum numeric peak2gene FDR to return
)
as_data_frame(p2g) %>% arrange(desc(Correlation)) %>%
head %>% knitr::kable(caption = "Peak-to-gene linkage")
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
| idxATAC | idxRNA | Correlation | FDR | VarQATAC | VarQRNA |
|---|---|---|---|---|---|
| 76211 | 2257 | 0.9776989 | 0 | 0.8779933 | 0.6813612 |
| 42485 | 14306 | 0.9621724 | 0 | 0.9966034 | 0.9910757 |
| 23457 | 10826 | 0.9599931 | 0 | 0.9372381 | 0.5969572 |
| 81836 | 3162 | 0.9599560 | 0 | 0.9999645 | 0.9820440 |
| 30183 | 12832 | 0.9559511 | 0 | 0.2635277 | 0.5640019 |
| 7959 | 1140 | 0.9556089 | 0 | 0.9458961 | 0.3378851 |
In the metadata the peakSet with ranges, strand and chromosome are saved, as well as the gene set with ranges, strand, chromosome and gene name.
metadata(p2g)
## $peakSet
## GRanges object with 141025 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 752506-753006 *
## [2] chr1 762688-763188 *
## [3] chr1 773653-774153 *
## [4] chr1 801007-801507 *
## [5] chr1 805041-805541 *
## ... ... ... ...
## [141021] chrX 154840753-154841253 *
## [141022] chrX 154842390-154842890 *
## [141023] chrX 154862033-154862533 *
## [141024] chrX 154912392-154912892 *
## [141025] chrX 154997007-154997507 *
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
##
## $geneSet
## GRanges object with 18601 ranges and 2 metadata columns:
## seqnames ranges strand | name idx
## <Rle> <IRanges> <Rle> | <array> <array>
## [1] chr1 69091 * | OR4F5 1
## [2] chr1 762902 * | LINC00115 2
## [3] chr1 812182 * | FAM41C 3
## [4] chr1 860530 * | SAMD11 4
## [5] chr1 894679 * | NOC2L 5
## ... ... ... ... . ... ...
## [18597] chrX 154299695 * | BRCC3 708
## [18598] chrX 154444701 * | VBP1 709
## [18599] chrX 154493852 * | RAB39B 710
## [18600] chrX 154563986 * | CLIC2 711
## [18601] chrX 154842622 * | TMLHE 712
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
##
## $seATAC
## [1] "/omics/groups/OE0533/internal/katharina/scDoRI/ArchR/Save_Vignette_object_4/Peak2GeneLinks/seATAC-Group-KNN.rds"
##
## $seRNA
## [1] "/omics/groups/OE0533/internal/katharina/scDoRI/ArchR/Save_Vignette_object_4/Peak2GeneLinks/seRNA-Group-KNN.rds"
atac_knn <- readRDS("Save_Vignette_object_5/Peak2GeneLinks/seATAC-Group-KNN.rds")
rna_knn <- readRDS("Save_Vignette_object_5/Peak2GeneLinks/seRNA-Group-KNN.rds")
atac_knn
## class: RangedSummarizedExperiment
## dim: 141025 491
## metadata(1): KNNList
## assays(2): ATAC RawATAC
## rownames(141025): f1 f2 ... f141024 f141025
## rowData names(13): score replicateScoreQuantile ... idx N
## colnames: NULL
## colData names(0):
dim(assays(atac_knn)[[1]])
## [1] 141025 491
#assays(atac_knn)[[1]] %>% head
#assays(rna_knn)[[1]] %>% head
We can also return loops:
p2g_loop <- getPeak2GeneLinks(
ArchRProj = proj,
corCutOff = 0.45,
resolution = 1000, # decrease resolution for plotting
returnLoops = TRUE
)
p2g_loop[[1]] %>% head
## GRanges object with 6 ranges and 2 metadata columns:
## seqnames ranges strand | value FDR
## <Rle> <IRanges> <Rle> | <numeric> <numeric>
## [1] chr1 762500-812500 * | 0.450391 7.66477e-25
## [2] chr1 762500-895500 * | 0.460307 4.97788e-26
## [3] chr1 762500-948500 * | 0.515106 2.34837e-33
## [4] chr1 812500-894500 * | 0.451671 5.41109e-25
## [5] chr1 812500-895500 * | 0.496871 9.11031e-31
## [6] chr1 812500-1004500 * | 0.460707 4.44943e-26
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
markerGenes <- c(
"CD34", #Early Progenitor
"GATA1", #Erythroid
"PAX5", "MS4A1", #B-Cell Trajectory
"CD14", #Monocytes
"CD3D", "CD8A", "TBX21", "IL7R" #TCells
)
p <- plotBrowserTrack(
ArchRProj = proj,
groupBy = "Clusters2",
geneSymbol = markerGenes,
upstream = 50000,
downstream = 50000,
loops = getPeak2GeneLinks(proj)
)
## GRanges object with 9 ranges and 2 metadata columns:
## seqnames ranges strand | gene_id symbol
## <Rle> <IRanges> <Rle> | <character> <character>
## [1] chr1 208059883-208084683 - | 947 CD34
## [2] chrX 48644982-48652717 + | 2623 GATA1
## [3] chr9 36838531-37034476 - | 5079 PAX5
## [4] chr11 60223282-60238225 + | 931 MS4A1
## [5] chr5 140011313-140013286 - | 929 CD14
## [6] chr11 118209789-118213459 - | 915 CD3D
## [7] chr2 87011728-87035519 - | 925 CD8A
## [8] chr17 45810610-45823485 + | 30009 TBX21
## [9] chr5 35856977-35879705 + | 3575 IL7R
## -------
## seqinfo: 24 sequences from hg19 genome
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
grid::grid.newpage()
grid::grid.draw(p$CD14)
The rows are clustered using k-means clustering based on value passed to parameter k (default = 25)
plotPeak2GeneHeatmap(proj, groupBy = "Clusters2")
## Warning: did not converge in 10 iterations